/*
 * FAtoms.cpp
 *
 *  Created on: May 13, 2012
 *      Author: marchi
 */

#include "FAtoms.h"
class brot{
	matrix co;
public:

	const matrix & operator()(const double a,const double b,const double c,
			const double alfa,const double beta,const double gamma){
	double degrad=M_PI/180.0;

	double ax=a;
    double alf=cos(degrad*alfa);
    double bet=cos(degrad*beta);
    double qt=sin(degrad*gamma);
    double gam=cos(degrad*gamma);
    double bx=b*gam;
    double by=b*qt;
    double cx=c*bet;
    double cy=c*(alf-bet*gam)/qt;
    double cz=sqrt(c*c-cx*cx-cy*cy);
    co[YY][XX]=0.0;
    co[ZZ][XX]=0.0;
    co[ZZ][YY]=0.0;
    co[XX][XX]=ax;
    co[XX][YY]=bx;
    co[XX][ZZ]=cx;
    co[YY][YY]=by;
    co[YY][ZZ]=cy;
    co[ZZ][ZZ]=cz;
    return co;
	}
};
void FAtoms::ReadaStep(FstreamC * fin){
	t_fileio *xd = fin->getfin();
	rvec   *x;
	matrix box;
	int    natos,step;
	real   prec,time;
	gmx_bool   bOK;
	int my_nframe=fin->gFrame();
	try{
		if(firsttime){
			if(xtc_seek_frame(fin->getfin(), my_nframe, nr)) throw string(" Beginning frame is out of range ");
			firsttime=false;
		}

		natos=nr;
		x=new rvec[nr];
		if(!read_next_xtc(xd,natos,&step,&time,box,x,&prec,&bOK))
			throw string("Cannot read next frame!");
	}
	catch(const string & s){
		cout << s << endl;
		exit(-1);
	}
	Metric Met(box);
	setCoord(Met,x);
	doCOtoOC();

	fin->nextFrame();
	delete [] x;
}

void FAtoms::ReadaStep(FstreamF * fin){
	ReadaStep(fin->getfin());
}

void FAtoms::ReadaStep(std::ifstream & fin){
	double table[6];
	fin.seekg(FORTRANBYTES,ios::cur);
	fin.read(reinterpret_cast<char *> (&table),sizeof(table));
	fin.seekg(FORTRANBYTES,ios::cur);
	table[0]*=unit_nm;
	table[2]*=unit_nm;
	table[5]*=unit_nm;
	double alpha = 90.0 - asin(table[1]) * 90.0 / M_PI_2;
	double beta = 90.0 - asin(table[3]) * 90.0 / M_PI_2;
	double gamma = 90.0 - asin(table[4]) * 90.0 / M_PI_2;

	Metric Met(brot()(table[0],table[2],table[5],alpha,beta,gamma));
	vector<float> X[DIM];

	for(int i=0;i<DIM;i++){
		X[i].resize(nr,0.0);
		fin.seekg(FORTRANBYTES,ios::cur);
		fin.read(reinterpret_cast<char *> (&X[i][0]),sizeof(X[0][0])*nr);
		fin.seekg(FORTRANBYTES,ios::cur);
	}
	rvec * x=new rvec[nr];
	for(int n=0;n<nr;n++){
		for(int o=0;o<DIM;o++) x[n][o]=X[o][n]*unit_nm;
	}

	setCoord(Met,x);
	doCOtoOC();
}
void FAtoms::moveOffset(FstreamC * fin){
	fin->nextFrame();
}
void FAtoms::moveOffset(FstreamF * fin){
	moveOffset(fin->getfin());
}
void FAtoms::moveOffset(std::ifstream & fin){
	const int DOUBLE=8;
	const int FLOAT=4;
	ios::streamoff OFFSET=FORTRANBYTES*2+6*DOUBLE+(FORTRANBYTES*2+nr*FLOAT)*DIM;
	fin.seekg(OFFSET,ios::cur);
}
FAtoms::~FAtoms() {
	// TODO Auto-generated destructor stub
}

